Timeseries clustering

Time series clustering is to partition time series data into groups based on similarity or distance, so that time series in the same cluster are similar.

Methodology followed:

• Use Variational Recurrent AutoEncoder (VRAE) for dimensionality reduction of the timeseries
• To visualize the clusters, PCA and t-sne are used

Paper: https://arxiv.org/pdf/1412.6581.pdf

Contents

1. Functions
2. Import
3. Load data and preprocess
4. Initialize VRAE object
5. Fit the model onto dataset
6. Transform the input timeseries to encoded latent vectors
7. Save the model to be fetched later

Functions

# https://gist.github.com/Zsailer/5d1f4e357c78409dd9a5a4e5c61be552
from IPython.display import HTML
from IPython.display import display

tag = HTML('''<script>
code_show=true; 
function code_toggle() {
    if (code_show){
        $('div.cell.code_cell.rendered.selected div.input').hide();
    } else {
        $('div.cell.code_cell.rendered.selected div.input').show();
    }
    code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
To show/hide this cell's raw code input, click <a href="javascript:code_toggle()">here</a>.''')
display(tag)


def load_one_dataset(direc = 'data', dataset="EMG", filename = "."):
    """
    Load one dataset.
    """
    datadir = direc + '/' + dataset + '/' + filename
    data = np.loadtxt(datadir, delimiter=',')
    data = np.concatenate((data, ), axis=0)
    data = np.expand_dims(data, -1)
    
    return data[:, 1:, :], data[:, 0, :]


def load_data(direc = 'data', dataset="EMG", all_file = [], do_pca = False, single_channel = None, 
             batch_size = 32, seq_len = 10, pca_component = 6):
    """
    Load all dataset and preprocess.
    
    all_file: list of all the filenames to load. Each have a segments * (seq length * channel amount) * 1 ndarray.
    do_pca: preform PCA on training data. If true, return original and transformed data, and PCA object.
            If False, return original data.
    single_channel: a list of channels to load, channel starts with 1. If None, load all channels.
    seq_length: should be the same as when converting to vrae format.
    
    Return:
    X_train: data in segments * seq length * features ndarray.
    y_train: labels in segments * 1 ndarray.
    X_train_ori: X_train before doing PCA.
    X_pca: PCA object.
    
    """
    
    # Load every dataset in list
    X_train = []
    y_train = []
    for file in all_file:
        X_train_small, y_train_small = load_one_dataset(direc = 'data', dataset="EMG", filename = file)
        X_train.append(X_train_small)
        y_train.append(y_train_small)
        print(f'Loading {file}, X shape {X_train_small.shape}, y shape {y_train_small.shape}', end = '')
        print(f', has label {np.unique(y_train_small)}')
    
    # Concatenate into one np array
    if len(all_file) == 1:
        X_train = X_train[0]
        y_train = y_train[0]
    else:
        X_train = np.concatenate(X_train, axis = 0)
        y_train = np.concatenate(y_train, axis = 0)
        
    # Cut the last several segments to match batch size
    X_train = X_train.reshape(X_train.shape[0], seq_len, -1)
    num_seg = (X_train.shape[0] // batch_size) * batch_size
    X_train = X_train[:num_seg, :, :]
    y_train = y_train[:num_seg, :]
    
    # Check
    if do_pca and single_channel:
        raise ValueError("Don't do both pca and single channel.")
    
    # Doing pca
    if do_pca:
        print('Doing PCA')
        # Copy original training data
        X_train_ori = np.copy(X_train)
        
        # Explained variance
        temp = X_train.reshape(-1, X_train.shape[2])
        X_pca = PCA(n_components=15).fit(temp)
        print(f'Explained variance ratio: {np.cumsum(X_pca.explained_variance_ratio_)}')
        
        # Need to specify n_components inside pca for later reconstruction
        X_pca = PCA(n_components = pca_component).fit(temp)
        X_train = X_pca.transform(temp)
        X_train = X_train.reshape(-1, seq_len, pca_component)
    
    # Extract single channels
    if single_channel:
        print(f'Extracting channels {single_channel}')
        single_channel = np.array(single_channel)-1
        X_train = X_train[:, :, single_channel]
        #X_train = np.expand_dims(X_train, -1)
    
    print(f'Dataset shape: {X_train.shape}')
    print(f'Label: {np.unique(y_train)}, shape: {y_train.shape}')
    
    if do_pca:
        return X_train, y_train, X_train_ori, X_pca
    
    return X_train, y_train


def recon(model, dataset):
    """
    Pass dataset through vrae to get a reconstruction.
    
    model: trained vrae model.
    dataset: original data in segments * seq length * features ndarray.
    
    return:
    reconstruction: segments * seq length * features ndarray.
    """
    
    torch_data = TensorDataset(torch.from_numpy(dataset))
    reconstruction = vrae.reconstruct(torch_data)
    reconstruction = reconstruction.transpose((1,0,2))
    
    return reconstruction


def plot_recon_feature(dataset, reconstruction, idx = None):
    """
    Plot the original and reconstructed feature of one segment.
    
    dataset: original data in segments * seq length * features ndarray.
    reconstruction: reconstructed data in segments * seq length * features ndarray.
    idx: Index of segment to plot. If None, randomly choose one.
    """
    
    num_seq = dataset.shape[0]
    num_features = dataset.shape[2]
    num_rows = -(-num_features//5)
    
    # reconstruction = recon(model, dataset)
    
    if idx:
        idx = idx-1
    else:
        idx = np.random.choice(num_seq, 1)[0]-1
    
    fig, axs = plt.subplots(num_rows, 5, figsize = (20, num_rows*5))
    for ii in range(num_features):
        ori = dataset[idx, :, ii]
        rec = reconstruction[idx, :, ii]
        axs[ii//5, ii%5].plot(ori, color = 'black')
        axs[ii//5, ii%5].plot(rec, color = 'red')
        axs[ii//5, ii%5].set_title(f'Feature #{ii+1}')
    fig.suptitle(f'Ori and rec of sequesce # {idx+1}', size = 20)
    plt.show()

    
def plot_recon_metrics(dataset, reconstruction, x_lim = None, verbose = True):
    """
    Plot correlation, mse, mean.
    
    dataset: original data in segments * seq length * features ndarray.
    reconstruction: reconstructed data in the same shape as dataset.
    
    """
    
    if x_lim:
        dataset = dataset[x_lim[0]:x_lim[1], :]
        reconstruction = reconstruction[x_lim[0]:x_lim[1], :]
    
    num_seq = dataset.shape[0]
    num_features = dataset.shape[2]

    mse_all = ((dataset-reconstruction)**2).mean(axis=1)
        
    mean_all = np.mean(dataset, axis = 1)
    

    corr_all = []
    for ii in range(num_features):
        corr_channel = []
        for jj in range(num_seq):
            corr_seq = np.corrcoef(reconstruction[jj, :, ii], dataset[jj, :, ii])[0,1]
            corr_channel.append(corr_seq)
        corr_channel = np.array(corr_channel)
        corr_all.append(corr_channel)
    corr_all = np.array(corr_all).transpose()
    
    if num_features == 1:
        fig, axs = plt.subplots(2, 1, figsize = (20, 12))
    else:
        fig, axs = plt.subplots(num_features, 1, figsize = (20, num_features*6))
    
    times = np.max(mse_all, axis = 0)
    for ii in range(num_features):
        axs[ii].plot(corr_all[:,ii]*times[ii]/3, color = 'r', label = 'corr')
        axs[ii].plot(mse_all[:, ii], color = 'y', label = 'mse')
        axs[ii].plot(mean_all[:, ii], color = 'dimgray', label = 'mean', alpha = 0.7)
        axs[ii].set_title(f'# {ii+1}, mean corr = {np.mean(corr_all[:,ii]):.4f}, ' \
                          f'mean mse = {np.mean(mse_all[:, ii]):.4f}, ' \
                          f'mean = {np.mean(mean_all[:, ii]):.4f}')
#         if x_lim:
#             axs[ii].set_xlim(x_lim)
        axs[ii].legend()
    if verbose:
        corr_mean = np.mean(corr_all, axis = 0)
        mse_mean = np.mean(mse_all, axis = 0)
        mean_mean = np.mean(mean_all, axis = 0)
        
        for jj in range(num_features):
            print(f'Channel {jj+1}, corr = {corr_mean[jj]:.4f}, '\
                  f'mse = {mse_mean[jj]:4f}, mean = {mean_mean[jj]:.4f}.')

        
def pca_inverse(PCA_obj, reconstruction):
    """
    Convert reconstructed PCs back to channels.
    
    PCA_obj: fitted PCA.
    reconstruction: in segments * seq length * features ndarray.
    """
    
    seq_len = reconstruction.shape[1]
    num_features = reconstruction.shape[2]
    
    reconstruction = reconstruction.reshape(-1, num_features)
    recon_channel = PCA_obj.inverse_transform(reconstruction)
    recon_channel = recon_channel.reshape(-1, seq_len, 15)
    
    return recon_channel


def visualize(z_run, y, inv_bhvs, one_in = 4, perplexity=80, n_iter=3000):
    """
    Visualize latent space using PCA and tSNE
    
    z_run: latent values, n_segments * lateng length ndarray.
    y: label of each segment.
    inv_bhvs: dict from label to behavior classes.
    one_in: default use one in every 4 segments.
    perplexity, n_iter: for tSNE.
    """
    
    z_run_down = z_run[::one_in, :]
    label = y[::one_in, :]
    
    # z_run_pca = TruncatedSVD(n_components=2).fit_transform(z_run_down)
    z_run_pca = PCA(n_components=2).fit_transform(z_run_down)
    z_run_tsne = TSNE(perplexity=perplexity, min_grad_norm=1E-12, n_iter=n_iter).fit_transform(z_run_down)
    
    all_colors = ['b','g','r','c','m','y','darkgrey']
    fig, axs = plt.subplots(1,2, figsize=(20,10))

    for ii in np.unique(label):
        ii = int(ii)

        x_pca = z_run_pca[:,0].reshape(-1,1)[label == ii]
        y_pca = z_run_pca[:,1].reshape(-1,1)[label == ii]
        axs[0].scatter(x_pca, y_pca, c=all_colors[ii], marker='.', label = inv_bhvs[ii], linewidths=None)

        x_tsne = z_run_tsne[:,0].reshape(-1,1)[label == ii]
        y_tsne = z_run_tsne[:,1].reshape(-1,1)[label == ii]
        axs[1].scatter(x_tsne, y_tsne, c=all_colors[ii], marker='.', label = inv_bhvs[ii], linewidths=None)

    axs[0].set_title('PCA on z_run')
    axs[1].set_title('tSNE on z_run')
    axs[0].legend()
    axs[1].legend()
    plt.show()

To show/hide this cell's raw code input, click here.

load_one_dataset \ load_data \ recon \ plot_recon_feature \ plot_recon_metrics \ pca_inverse \ visualize

Import required modules

from vrae.vrae import VRAE
from vrae.utils import *
import numpy as np
import torch
import pickle

from sklearn.decomposition import PCA, TruncatedSVD
from sklearn.manifold import TSNE
from sklearn.metrics import mean_squared_error as mse

import plotly
from torch.utils.data import DataLoader, TensorDataset
plotly.offline.init_notebook_mode()

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Download dir

dload = './model_dir'

Hyper parameters

seq_len = 10
hidden_size = 256
hidden_layer_depth = 3
latent_length = 32
batch_size = 32
learning_rate = 0.00002
n_epochs = 2000
dropout_rate = 0.0
optimizer = 'Adam' # options: ADAM, SGD
cuda = True # options: True, False
print_every=10
clip = True # options: True, False
max_grad_norm=5
loss = 'MSELoss' # options: SmoothL1Loss, MSELoss
block = 'LSTM' # options: LSTM, GRU
output = False

Load data and preprocess

training_file = ['20201020_Pop_Cage_001','20201020_Pop_Cage_002','20201020_Pop_Cage_003','20201020_Pop_Cage_004',
                '20201020_Pop_Cage_006']
X_train, y_train, X_train_ori, X_pca = load_data(direc = 'data', dataset="EMG", all_file = training_file,
                         do_pca = True, single_channel = None, 
                         batch_size = batch_size, seq_len = seq_len, pca_component = 6)
train_dataset = TensorDataset(torch.from_numpy(X_train))

Loading 20201020_Pop_Cage_001, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.]
Loading 20201020_Pop_Cage_002, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.]
Loading 20201020_Pop_Cage_003, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.  4.]
Loading 20201020_Pop_Cage_004, X shape (3601, 150, 1), y shape (3601, 1), has label [-1.  0.  1.  2.  3.  4.]
Loading 20201020_Pop_Cage_006, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.  4.]
Doing PCA
Explained variance ratio: [0.7165347  0.82274211 0.88277787 0.92401944 0.9436966  0.95823764
 0.97067536 0.9778943  0.98333348 0.98764237 0.99133662 0.99435286
 0.99654689 0.99843935 1.        ]
Dataset shape: (17984, 10, 6)
Label: [-1.  0.  1.  2.  3.  4.], shape: (17984, 1)

num_features = X_train.shape[2]

Initialize VRAE object

VRAE inherits from sklearn.base.BaseEstimator and overrides fit, transform and fit_transform functions, similar to sklearn modules

vrae = VRAE(sequence_length=seq_len,
            number_of_features = num_features,
            hidden_size = hidden_size, 
            hidden_layer_depth = hidden_layer_depth,
            latent_length = latent_length,
            batch_size = batch_size,
            learning_rate = learning_rate,
            n_epochs = n_epochs,
            dropout_rate = dropout_rate,
            optimizer = optimizer, 
            cuda = cuda,
            print_every=print_every, 
            clip=clip, 
            max_grad_norm=max_grad_norm,
            loss = loss,
            block = block,
            dload = dload,
            output = output)

/home/roton2/miniconda3/envs/emg/lib/python3.9/site-packages/torch/nn/_reduction.py:42: UserWarning:

size_average and reduce args will be deprecated, please use reduction='sum' instead.

Fit the model onto dataset

#vrae.fit(train_dataset)

#If the model has to be saved, with the learnt parameters use:
vrae.fit(train_dataset)

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
/tmp/ipykernel_3388461/2701181639.py in <module>
      2 
      3 #If the model has to be saved, with the learnt parameters use:
----> 4 vrae.fit(train_dataset)

~/ruize/timeseries-clustering-vae/vrae/vrae.py in fit(self, dataset, save)
    369             #print('Epoch: %s' % i)
    370 
--> 371             average_loss = self._train(train_loader)
    372 
    373             if (i + 1) % self.print_every == 0:

~/ruize/timeseries-clustering-vae/vrae/vrae.py in _train(self, train_loader)
    330             self.optimizer.zero_grad()
    331             loss, recon_loss, kl_loss, _ = self.compute_loss(X)
--> 332             loss.backward()
    333 
    334             if self.clip:

~/miniconda3/envs/emg/lib/python3.9/site-packages/torch/_tensor.py in backward(self, gradient, retain_graph, create_graph, inputs)
    253                 create_graph=create_graph,
    254                 inputs=inputs)
--> 255         torch.autograd.backward(self, gradient, retain_graph, create_graph, inputs=inputs)
    256 
    257     def register_hook(self, hook):

~/miniconda3/envs/emg/lib/python3.9/site-packages/torch/autograd/__init__.py in backward(tensors, grad_tensors, retain_graph, create_graph, grad_variables, inputs)
    145         retain_graph = create_graph
    146 
--> 147     Variable._execution_engine.run_backward(
    148         tensors, grad_tensors_, retain_graph, create_graph, inputs,
    149         allow_unreachable=True, accumulate_grad=True)  # allow_unreachable flag

KeyboardInterrupt: 

Plot loss and MSE

plt.plot(vrae.all_loss)

[<matplotlib.lines.Line2D at 0x7f140443bd60>]

plt.plot(vrae.rec_mse)

[<matplotlib.lines.Line2D at 0x7f1293797cd0>]

Transform the input timeseries to encoded latent vectors

#If the latent vectors have to be saved, pass the parameter `save`
z_run = vrae.transform(train_dataset, save = True, filename = 'z_run_e57pca_2000epoch.pkl')
z_run.shape

(17984, 32)

Save / load the model

vrae.save('./vrae_e5_3000epoch.pth')

vrae.load(dload+'/vrae_e57pca_2000epoch.pth')
with open(dload+'/z_run_e57pca_2000epoch.pkl', 'rb') as fh:
    z_run = pickle.load(fh)

Reconstruction

reconstruction = recon(vrae, X_train)

plot_recon_feature(X_train, reconstruction, idx = None)

plot_recon_metrics(X_train, reconstruction, x_lim = [2000, 4000])

Channel 1, corr = 0.9864, mse = 5.760487, mean = -1.2134.
Channel 2, corr = 0.9550, mse = 6.425530, mean = 0.7497.
Channel 3, corr = 0.9196, mse = 9.416685, mean = -4.8806.
Channel 4, corr = 0.9273, mse = 9.425540, mean = 3.9921.
Channel 5, corr = 0.8631, mse = 10.365440, mean = 1.4297.
Channel 6, corr = 0.8351, mse = 13.525274, mean = 3.2135.

Convert PC back to channel traces

recon_channel = pca_inverse(X_pca, reconstruction)

plot_recon_feature(X_train_ori, recon_channel, idx = None)

plot_recon_metrics(X_train_ori, recon_channel, x_lim = [0, 2000])

Channel 1, corr = 0.7173, mse = 46.461902, mean = 25.5160.
Channel 2, corr = 0.6914, mse = 41.082651, mean = 24.5046.
Channel 3, corr = 0.6536, mse = 62.893647, mean = 28.5684.
Channel 4, corr = 0.5966, mse = 43.409875, mean = 19.7128.
Channel 5, corr = 0.6227, mse = 23.129086, mean = 12.4884.
Channel 6, corr = 0.7438, mse = 11.456219, mean = 30.5044.
Channel 7, corr = 0.8457, mse = 43.196223, mean = 46.2902.
Channel 8, corr = 0.8658, mse = 49.333910, mean = 49.2048.
Channel 9, corr = 0.6768, mse = 27.816960, mean = 18.7294.
Channel 10, corr = 0.7323, mse = 26.337528, mean = 24.9044.
Channel 11, corr = 0.8376, mse = 4.954293, mean = 39.9276.
Channel 12, corr = 0.5631, mse = 87.285175, mean = 17.8411.
Channel 13, corr = 0.8585, mse = 76.566824, mean = 44.3026.
Channel 14, corr = 0.8356, mse = 64.866368, mean = 38.6561.
Channel 15, corr = 0.7998, mse = 48.271006, mean = 34.4277.

Reconstruction of test dataset

testing_file = ['20201020_Pop_Cage_005', '20201020_Pop_Cage_007']
X_test, y_test, X_test_ori, test_pca = load_data(direc = 'data', dataset="EMG", all_file = testing_file,
                         do_pca = True, single_channel = None, 
                         batch_size = batch_size, seq_len = seq_len, pca_component = 6)

Loading 20201020_Pop_Cage_005, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.  5.]
Loading 20201020_Pop_Cage_007, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.  4.]
Doing PCA
Explained variance ratio: [0.72963202 0.82936616 0.88681766 0.92899866 0.9469495  0.9610898
 0.97256296 0.97950278 0.98456741 0.98865807 0.99216837 0.99489398
 0.99686689 0.99867192 1.        ]
Dataset shape: (7168, 10, 6)
Label: [-1.  0.  1.  2.  3.  4.  5.], shape: (7168, 1)

recon_test = recon(vrae, X_test)
recon_channel_test = pca_inverse(test_pca, recon_test)

plot_recon_feature(X_test, recon_test, idx = None)

plot_recon_metrics(X_test, recon_test, x_lim = [0, 2000], verbose = True)

Channel 1, corr = 0.9468, mse = 225.725982, mean = 3.9363.
Channel 2, corr = 0.8644, mse = 183.498191, mean = 0.3208.
Channel 3, corr = 0.8020, mse = 106.458381, mean = -0.2006.
Channel 4, corr = 0.8002, mse = 119.267475, mean = 1.7936.
Channel 5, corr = 0.7036, mse = 99.302438, mean = 0.7821.
Channel 6, corr = 0.6646, mse = 78.343952, mean = 1.4382.

recon_test_channel = pca_inverse(test_pca, recon_test)

plot_recon_feature(X_test_ori, recon_test_channel, idx = None)

plot_recon_metrics(X_test_ori, recon_test_channel, x_lim = [0, 2000])

Channel 1, corr = 0.7253, mse = 95.631335, mean = 32.4938.
Channel 2, corr = 0.7342, mse = 53.446205, mean = 28.8226.
Channel 3, corr = 0.6596, mse = 99.875063, mean = 33.7234.
Channel 4, corr = 0.6062, mse = 47.280105, mean = 20.6514.
Channel 5, corr = 0.5979, mse = 34.336671, mean = 13.4626.
Channel 6, corr = 0.6655, mse = 119.873159, mean = 31.9634.
Channel 7, corr = 0.8184, mse = 134.456782, mean = 51.8214.
Channel 8, corr = 0.8430, mse = 130.019252, mean = 57.4305.
Channel 9, corr = 0.6882, mse = 44.186982, mean = 22.3521.
Channel 10, corr = 0.7299, mse = 76.777045, mean = 32.1617.
Channel 11, corr = 0.7891, mse = 184.603463, mean = 47.7612.
Channel 12, corr = 0.5422, mse = 145.114475, mean = 21.8368.
Channel 13, corr = 0.8438, mse = 171.366514, mean = 53.2006.
Channel 14, corr = 0.7937, mse = 120.049905, mean = 42.0137.
Channel 15, corr = 0.7714, mse = 115.757965, mean = 39.3028.

This is another model without pca

testing_file = ['20201020_Pop_Cage_005']
X_test, y_test = load_data(direc = 'data', dataset="EMG", all_file = testing_file,
                         do_pca = False, single_channel = None, 
                         batch_size = batch_size, seq_len = seq_len, pca_component = 6)

Loading 20201020_Pop_Cage_005, X shape (3599, 150, 1), y shape (3599, 1), has label [-1.  0.  1.  2.  3.  5.]
Dataset shape: (3584, 10, 15)
Label: [-1.  0.  1.  2.  3.  5.], shape: (3584, 1)

recon_test = recon(vrae, X_test)

plot_recon_feature(X_test, recon_test, idx = None)

plot_recon_metrics(X_test, recon_test, x_lim = [0, 2000])

Visualize latent space using PCA and tSNE

bhvs = {'crawling': np.array([0]),
        'high picking treats': np.array([1]),
        'low picking treats': np.array([2]),
        'pg': np.array([3]),
        'sitting still': np.array([4]),
        'grooming': np.array([5]),
        'no_behavior': np.array([-1])}

inv_bhvs = {int(v): k for k, v in bhvs.items()}

visualize(z_run, y = y_train, inv_bhvs = inv_bhvs)

test_dataset = TensorDataset(torch.from_numpy(X_test))
z_run_test = vrae.transform(test_dataset, save = False)
visualize(z_run_test, y = y_test, inv_bhvs = inv_bhvs, one_in = 4)

# # Create clusters.annot
# from sklearn.mixture import GaussianMixture

# # Predict cluster assignments
# gm = GaussianMixture(n_components=5, random_state=0).fit(z_run)
# clusters = gm.predict(z_run)

# # Number of seconds in each sequence
# filt_time_step = 0.025
# num_secs_seq = sequence_length * filt_time_step
# end_time = len(z_run) * num_secs_seq + num_secs_seq

# # Print head of the file
# f = open ('Pop01-06_18_2021.annot','w')
# # write the header--------------------
# f.write('Bento annotation file\n')
# f.write('Movie file(s): {}\n\n'.format('Pop_20210618_cage_C1_01.avi'))
# f.write('{0} {1}\n'.format('Stimulus name:',''))
# f.write('{0} {1}\n'.format('Annotation start frame:',1))
# f.write('{0} {1}\n'.format('Annotation stop frame:', 26994))
# f.write('{0} {1}\n'.format('Annotation framerate:', 30))

# f.write('\n{0}\n'.format('List of channels:'))
# channels = ['cluster_num']
# for item in channels:
#         f.write('{0}\n'.format(item))
# f.write('\n');

# f.write('{0}\n'.format('List of annotations:'))
# clust_names = ['cluster_{}'.format(str(num)) for num in set(clusters)]
# labels = clust_names
# # labels = [item.replace(' ','_') for item in labels]
# for item in labels:
#     f.write('{0}\n'.format(item))
# f.write('\n')

# # now write the contents---------------
# for ch in channels:
#     f.write('{0}----------\n'.format(ch))
#     for beh in labels:
#         f.write('>{0}\n'.format(beh))
#         f.write('{0}\t {1}\t {2} \n'.format('Start','Stop','Duration'))

#         idxs = np.where(clusters == int(beh.split('_')[-1]))[0]
#         for hit in idxs:
#             start_time = hit * num_secs_seq/2
#             end_time = start_time + num_secs_seq
#             f.write('{0}\t{1}\t{2}\n'.format(start_time, end_time, num_secs_seq))
#         f.write('\n')
#     f.write('\n')

# f.close()